Measuring Renyi Entanglement Entropy with Quantum Monte Carlo 
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We develop a quantum Monte Carlo procedure, in the valence bond basis, to measure the Renyi 
entanglement entropy of a many-body ground state as the expectation value of a unitary Swap 
operator acting on two copies of the system. An improved estimator involving the ratio of Swap 
operators for different subregions enables convergence of the entropy in a simulation time polyno- 
mial in the system size. We demonstrate convergence of the Renyi entropy to exact results for a 
Heisenberg chain. Finally, we calculate the scaling of the Renyi entropy in the two-dimensional 
Heisenberg model and confirm that the Neel groundstate obeys the expected area law for systems 
up to linear size L = 32. 
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The measurement of entanglement entropy provides 
new tools to determine universal properties of interact- 
ing quantum many-body systems in condensed-matter 
physics. For example, conformally invariant systems in 
one dimension (ID) display universal logarithmic scaling 
of entanglement entropy [ij. Many quantum systems in 
two dimensions (2D) and higher are predicted to have 
"area law" scaling of entanglement entropy, with either 
universal additive logarithmic corrections for critical sys- 
tems [3, S] or universal additive constant corrections for 
topologically ordered systems 0, Q . 

Scalable simulation methods have become crucial for 
the study of groundstate, finite-temperature, and criti- 
cal properties of quantum many-body systems. Unfortu- 
nately, there are no known scalable simulation methods 
to calculate entanglement entropy in D > 1. Calculating 
the von Neumann entropy (S'l, defined below) requires 
calculating the reduced density matrix, which can be 
done in Lanczos or density matrix renormalization group 
(DMRG). However, scalable simulation techniques in di- 
mensions D > I, are so far essentially restricted to quan- 
tum Monte Carlo (QMC), which does not provide access 
to the reduced density matrix, but only approximates it 
via a Metropolis importance-sampling algorithm. 

We resolve this problem and show that a generalized 
Renyi entropy (5*2, defined below) can be related to the 
expectation value of a unitary Swap operator in the va- 
lence bond basis, directly accessible in QMC. Using QMC 
simulations, we calculate iS'2 for a spin 1/2 Heisenberg 
model on ID and 2D lattices. In ID, we demonstrate 
that ^2 calculated with QMC converges to the exact re- 
sult obtained with DMRG simulations. In 2D, we show 
that the leading-order term in the scaling of S'2 goes like 
the boundary of a subregion A, confirming the area law 
expected from recent studies [6j. 

The generalized Renyi entanglement entropies are de- 
fined by 




Snip/ 



■In [Tr(pl)], 



(1) 



SwapA\Vr] 



FIG. 1: (color online) A six-site chain, with two non- 
interacting copies (top and bottom) before (a) and after (b) 
the SwapA operation. The region A consists of three light- 
colored sites on the left; the complement region B of three 
dark sites on the right. The curved lines denote singlets in 
the state \ Vr}, which is a product of two different valence bond 
states, one per copy. The ground state of the entire system is 
a linear combination of similar \Vr)- 



where pA is the reduced-density matrix of a subregion 
A entangled with its complement B, and n > 0. In the 
limit n — ^ 1, one recovers the familiar von Neumann en- 
tropy, Si{pa) = — Tr(p^ \npA)- Recently, the generalized 
Renyi entropies have attracted considerable attention in 
the condensed matter community, due to their ability to 
encode information about the whole "entanglement spec- 
trum" of pA^ together allowing the set of Sn to contain 
much more information than S'l alone 0- For example, 
the concept of topological entanglement entropy has re- 
cently been generalized to the family of Renyi entropies, 
where it has been shown to be equal to the logarithm of 
the total quantum dimension, S'top oc log!?, independent 
of n Q. Universal corrections to the scaling of Sn have 
also been calculated in a field-theoretic treatment of the 
0{N) model Q. Further, any two Renyi entropies 5'„ 
and Sm obey the inequality Sn > Sm for n < m, making 
S'2 a useful lower bound on Si. 

Valence bond basis and the Swap operator - Wc briefly 
review the notation of Sandvik, and refer the reader to 
Refs. [IMll for additional details on the valence bond 
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basis. One begins by writing the expansion of a singlet 
wavefunction of N spins as N/2 valence bonds: 

\^o) = E MK'bl) . . . 65v/2)> = E MVr) (2) 

r r 

where the valence bonds are singlet pairs denoted by 
(a, 6) = (Itaib) — \iatb))/V^, occurring between sites a 
and b of the opposite sublattices of a bipartite lattice, 
and r labels all possible tilings of bonds. The coefficients 
fr are unknown, but an importance sampling scheme [lo| 
for the valence bond state \Vr) will be outlined below. 

To derive an estimator for 5*2 , we construct two copies 
of the system (higher entropies Sn require n copies) ; one 
"real", one "ancillary" (Fig.[T]). Then we define an oper- 
ator Swap A which acts between the copies of the system, 
swapping the configurations within region A. To define 
this operator, we give its matrix elements in a product 
basis, e.g. the 5^ basis. Let |a) be a complete basis 
of states in the region A and let be a complete ba- 
sis of states in the complement region B. The state of 
each copy can be decomposed in this product basis as 
I*) = /3 C'a.^lof)!/?) for some amplitudes Ca,p- We 
define the unitary operator SwapA by 

SwapA^E C'ai,/3ll"l)l/3l)) ® (E ■Da2,/?2l"2)|/32)) 

Ql.,3l a2,P2 

= E ^"lA E ^"2.&(l«2>|/3i>) ® (|ai)|/32)). (3) 

Ql,/3l 012,02 

The two copies of the system will not interact with each 
other, so the ground state of the combined system will be 
the product of the ground states on each copy: 1 4'o CS" ^'o) ■ 
Thus, the expectation value of SwapA is 

= Cai,l3iCa2,PiCa2, 02^01,02 

Ql,Q2,/3l,^2 

= E (^^)"i'"2(PA)Q2,ai Tr(p^), (4) 

ai ,a2 

where (pa)qi.q2 = J^fSi Cai,i3iC a2,Pi denotes a matrix 
element of pA ■ Therefore, 

S2{pa) = - ln(Tr(p^)) = - lii{{SwapA))- (5) 

Equation ([5]) is basis independent In particular it 

holds in the valence bond basis, where the SwapA oper- 
ator is not defined by Eq. dSj. Rather, it is defined to 
swap the endpoints of valence bonds between the real and 
ancillary copies in the region A, as illustrated in Fig. [T] 
Measuring the Swap operator in QMC- We now 
present a procedure for measuring the entropy 5*2 via 
importance sampling of the Swap a operator in QMC. A 
projector scheme in the valence bond basis has recently 
been pioneered by Sandvik, and we only briefly review 



the notation here, referring the reader to Refs. 
for details of implementation. The method is a T = 
projector QMC, where high powers of a Hamiltonian H 
are used to project out the ground state |^o) from a trial 
wavefunction ^': cx |*o)- In Sandvik's QMC 

scheme, (— i/)™ is written as a sum of all products of m 
bond operators, 

m 

(-i/r ^EH^^T^E^- (6) 

r 1=1 r 

where for concreteness we use the spin 1/2 Hcisenbcrg 
model, defining H = - J2{a.b) ^^b and Hat = -(Sa • S;, - 
1/4). The "operator string" is sampled according to 
its weight, Wr, accrued upon evolution of a trial valence 
bond state \V) under projection: 

Pr\V) = Wr\V{r)). (7) 

As shown in Ref. for the Hcisenbcrg model, Wr is 
simply related to the number of off-diagonal operations 
77ioff in the projection, Wr = 2^™°'*'. 

To sample the Swap a operator, one req uires a double- 
projector valence bond QMC scheme [lO|, where a gen- 
eral expectation value for any observable O is given by 

^ Y.rl{y\P:OPr\V) ^ Y.rlWlWr{V {l)\0\V {r)) 
^ ^ T.rmPtPr\V) i:rlWlWr{V{l)\V{r)) ' 

(8) 

In this case, the two operator strings. Pi and P^, are ap- 
plied to two copies of the system; the expectation value 
of the Swap A operator as illustrated in Fig. [1] can then 
be calculated directly. Specifically, one performs im- 
portance sampling of operator strings according to the 
weight WiWr{V {l)\V [r)) , and measures the QMC aver- 
age expectation value 

/ {V{l)\SwapA\V{r)) \ 

^^^-P^) = \ (y{i)\v{r)) /' 

calculating S2 from Eq. ([5]). 

Results for a ID chain are illustrated in Fig. [2] There, 
simulations were performed using a double-projector 
QMC, with a simple "columnar" trial state \ V) (alternat- 
ing nearest-neighbor bonds in ID). In the following data, 
four random bond operators were changed in each oper- 
ator string per Monte Carlo step, and (unless otherwise 
stated) m = 407V. Exact results for S2 were obtained 
with a DMRG simulation, converged with 1000 states. 
One immediately sees that the naive expectation value 
(Swopa) results in very large statistical error bars when 
the region A grows large (in ID, A is the linear region of 
size i). One may understand this by considering the ex- 
pected scahng of 6*2. Namely, 5*2 scales logarithmically 
with the size of A, so the expectation value (SwapA) 
should be polynomially small in A. As a result, there are 
only very few configurations contributing significantly to 
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FIG. 2: (color online) The Renyi entropy S2 as a function of 
site index i £ A, for a 100-site Heisenberg chain with open 
boundaries, calculated with DMRG and QMC. Data labeled 
"Swap" was calculated with Eq. ([9]) with one QMC simula- 
tion, while data labeled jmax = 5 was calculated with Eq. (|10|) 
using 20 separate QMC simulations with a range of j G [1,5]. 
The inset shows the convergence of S2 to the exact value 
(dashed line) for i — 6 with up to m = 4000. 



{Swap a), so in practice the importance sampling is done 
poorly, resulting in large error bars and possibly jeopar- 
dizing simulation ergodicity [16[. The problem is even 
worse in 2D and higher, where the expectation value 
(SwapA) should be exponentially small in A. 

To combat this issue, we propose a refinement to the 
algorithm, which we called "improved ratio" sampling. 
Consider first a one-dimensional chain, and define regions 
A^ , A^ , A" , such that A^ contains i sites, and A'^ is 
a subset of In other words, each region is 

obtained by adding one site to region A^, and region A" 
is the empty set (thus, SwapAo is equal to the identity 
operator). In a given simulation, one can calculate the 
ratio 



jSwapA^^i) ^ Eri WiWr{V{l)\SwapA^^AV{'r)) 
{SwapA,) Eri WiWr{V{l)\SwapAAV{r)) 



(10) 



where, for each i = 0, ...,n— 1, the log of this ratio is equal 
to minus the difference S2{pa'+^) ~ S2{pa^)- Computing 
this ratio for each i will let us compute 6*2 for all A*. To 
calculate this ratio, the simulation must be performed 
with the modified sampling weight. 



WiWr{V{l)\SwapAAV{r)), 



(11) 



i.e., a unique QMC simulation must be done for each de- 
sired i . If every site were used for a different simulation 
weight, the QMC algorithm would obtain an additional 
multiplicative factor of N in its scaling. In practice, this 
can be reduced by noting that good statistical control 
can be retained by calculating (SwapAi+j) / (SwapAi) for 
fixed i and a range of j S [1, jmax] (see also the Discus- 
sion). We illustrate this in Fig. [51 where simulations with 
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FIG. 3: (color online) The Renyi entropy divided by the linear 
dimension £ of the entangled region A, for 2D lattices with 
L = 8 and L — 16. The data labeled "Swap" is from a 
single simulation calculating Eq. ((9| directly. The other data 
sets are derived from the improved ratio estimator, Eq. (|10[) . 
with different ranges of r G [l,rmax] (see text). The inset is 
a periodic L — 8 lattice with region A consisting of the 16 
central (dark) sites labeled by ^ = 4. 



jmax = 5 are sufficient to converge 100-site chain to the 
exact results along most of its length. 

2D Heisenberg Results - We now extend these concepts 
to QMC simulations of the spin 1/2 Heisenberg model on 
N = L X L lattices. For the simulations using the im- 
proved ratio estimator, {Swap A^+r-) / {Swap a') , we define 
A' as a square region of linear size I containing ^ x ^ sites, 
so A^'^'^ contains 2lr + more sites than A^. 

Figure [3] illustrates 5*2 on two 2D Heisenberg mod- 
els with L = 8 and L = 16. In the case of L = 8, 
data for a single simulation calculating {SwapA) directly 
is essentially identical to the improved ratio estimators 
for all I. Here, the direct estimator (Eq. ^) corre- 
sponds to one simulation, while rmax = 1 corresponds 
to 6 different simulations of the improved ratio estima- 
tor {Swap Ai+i) / {Swap At) and I e [1,6]. In contrast to 
L = 8, for L = 16 the direct estimator is significantly 
different than the improved ratio estimator. We can see 
the convergence of the data as one successively decreases 
the range of r G [l,rniax], so that rmax = 2 is identical 
to Tinax = 1, which is the smallest possible range for the 
improved ratio estimator in the 2D geometry illustrated. 
It is important to note that for large rmax the correct 
value of 5*2 does not lie within the error bars of the data. 
This is an indication that simulation ergodicity may be 
jeopardized by low sampling of the SwapA operator. 

We use the improved ratio estimator with rmax = 1 to 
scale the results for the Renyi entropy of the Heisenberg 
model to larger system sizes in 2D. Results are given in 
Fig. H for L = 4 to L = 32, plotted as S2/t From 
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FIG. 4: (color online) Scaling of the Renyi entropy divided by 
the linear dimension I of the entangled region A for different 
systems sizes in 2D. All data is calculated with the improved 
ratio estimator, Eq. (|10|) . with rmax = 1. 



Rcf. [6| one expects the scaling of S\ in the Neel state 
to obey the area law; also recall that S'2 < 5*1. The 
data illustrated is clearly consistent with the area law 
S^ji ^ const, for <C i, at which point boundary effects 
likely become important. In particular, multiplicative 
log corrections (which were apparent in similar system 
sizes and geometries for the Valence-Bond entanglement 
entropy 14, 1^) are not present. 

Discussion - In this paper we have presented an al- 
gorithm for measuring the Renyi entropy S'2 in valence 
bond basis QMC simulations via the expectation value of 
a swap operator between two copies of the system. Us- 
ing an improved ratio estimator, we are able to converge 
the expectation value of ^2 to the exact result on a ID 
Heisenberg chain. Using the same procedure, we have 
presented the first measurement of the Renyi entangle- 
ment entropy in a 2D system, confirming the area law for 
the Neel groundstatc of the Heisenberg model. 

The simple double-projector QMC algorithm used in 
this paper is known to scale approximately as O(to^), 
[1^ assuming m > N . Thus, the current simulation re- 
sults for the direct swap expectation value, Eq. (jH), also 
have 0{m?) scaling, while results obtained using the im- 
proved ratio estimator, Eq. PH)) . have 0{Lm?) scaling 
in the current geometries. We note however that these 
geometries may not be ideal for very large i and L, even 
in the improved ratio sampling. For example, for L = 32, 
results for rmax = 2 and rmax = 1, although remaining 
converged within error bars, begin to develop larger dis- 
crepancies when £ approaches L. 

From the current work (which used about 10 CPU- 
years), the question naturally arises whether it is possible 
to converge the expectation value of the swap operator on 
larger system sizes. We expect that this will be possible 
using a different sequence of regions in the improved 
ratio estimator. Rather than considering a sequence of 



regions A^ consisting of £-by-^ squares, we can use a more 
finely grained sequence of regions; e.g. by defining A^ to 
contain i sites such that for j = £^ , AMs an £-by-^ square. 
Incrementing the number of sites in i by one would result 
in (at worst) 0{Nrin?) scaling in the current algorithm. 
However, we also note that recently-developed loop algo- 
rithms result in significant improvement in scaling, allow- 
ing for the convergence of observables such as energy and 
spin correlation functions for L = 256 and larger [12l | . In 
these schemes, the scaling of the direct swap operator 
(Eq. dni)) would decrease to 0{m), while the expected 
scaling of the improved ratio estimator would be 0{Nm) 
at worst, depending on the definitions of the A"^ regions. 
We therefore expect that the algorithm will require only 
a polynomial number of samples to converge the Renyi 
entropy on arbitrary system sizes. 
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